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ABSTRACT 


Refraction arrivals from the upper mantle were 
mecerdedion 233500 km. profide stretching from, Lake 
Superior to Alaska. Recordings from the middle section 
of the profile (1000 - 2000 km.) were digitized, and 
various frequency and velocity filters were used to 
enhance this data. One version of the velocity filter, 

a non-linear operator, obtained the highest degree of 
data enhancement, but at the expense of signal distortion. 
The analysis presented in this paper rests primarily 

on these enhanced records and secondarily on the analog 


records, 


The velocity structure was obtained by matching 
a theoretical travel time curve calculated from proposed 
structures to the observed travel time curve. The most 
successful model thus produced was model U of A #3. 
This model contained four striking features; two velocity 
reversals, one at a depth of 42 to 70 km. and the other 
aul e> co. 70m, , Oth OL awuichewere nm ollowed by. raoid 
velocity increases, and two rapid velocity increases at 
depths of 455 and 660 km. These velocity reversals were 
necessary to reproduce the shadow zone at a distance of 
900 km, and its accompanying cusp at 300 km., and the 
shadow zone at 1500 km. with its cusp at 800 km. The 


rapid velocity increases were used to produce two cusp 
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bairs, »tie first avy 3200 and 1/00 km., and the second at 
3300 and 2200 km. Two additional velocity reversals 
aeecepuns Of 195 to 230 km, and 250 Lo 325 km. ath a 
small velocity increase between them were neede to 

produce the shadow zones at 2000 and 2400 km. respectively. 
These reversals were too small to produce any detectable 


Cus ping, 


TSPcomprece che Cetaiice Of Gils suructure, 
and to obtain the velocity behaviour in the transition 
zones, the theoretical amplitudes calculated from the 
model must be made to match the observed amplitudes. 
Thais amplitude study is not’ presented in this paper, 
However, the computer program listed in the appendix 
calculates both amplitudes and travel times, and is 


presently being used to complete this study. 
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DATA COLLECTION 


ane itncroduchi0n 


Project Early Rise, sponsored by the United 
staves Air Force Office of Scientific Research and co- 
ordinated by the United States Geological Survey, was 
aimeexperiment, designed fom the investigation of the 
upper mantle using seismic body waves. In previous 
studies, earthquakes and nuclear explosions were used as 
UheeScismucssounce, oweverm,. Cue slo ene nature sor .bhese 
sources, «the resulting refraction proiailes were often in 
anomalous areas and were necessarily short or sparcely 
occupied, Also, ache source, parameters, .such as ,origin 
bame, saypoOCcenterm, jenerey, and yuranster function, were snot 
always well known, which added further confusion to the 
complex travel time and amplitude curves for seismic body 
Waves. For a suceessful experiment,. the profiles had ito 
be long and detailed, with as many of the source parameters 
independently measurable or known as possible. Project 


Early Rise was planned to meet these requirements. 
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information been available for the planning of a refraction 
program. Cumming and Kanasewich (1965) demonstrated the 
feasibility of using underwater chemical explosions of a 
reasonable size as the seismic source for upper mantle 
investigation. They recorded several of the two to ten 
ton shots from the 1963 Lake Superior crustal experiment 

at distances sufficiently great (1600 - 1900 km.) that the 
seismic waves producing the first events had penetrated the 
mantle. This experiment also provided a known crust 

(Berry and West, 1966) under Lake Superior and provided 
data regarding seismic coupling of underwater explosions 

at various locations within the lake. Underwater shots had 
a second advantage due to the fact that all shots could be 
set off at the same location, thus yielding a common 
Cransther LUNnGLIOnN, coupling andrcrusy tOrvall explosions, 
The same experiment provided valuable experience in shot 
procedure. Most important, it resulted in a proven system 
of hydrophone stations for measuring the shot instant and 
precise shot location. The 1963 experiment, as well as 
another one in 1964, guaranteed a reliable source with 


accurately measurable source parameters. 


As a result of these studies, forty shots of five 
tons each were detonated at a common shot point in July 
1966 during Project Early Rise. These formed the sources 
for the many seismic profiles, shown in figure 1.1, under- 


taken by various university and government research 
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laboratories. Five seismic crews, two from the University 
of Alberta and three from the United States Geological 
Survey (USGS), employing mobile linear arrays, participated 
in the Western Canada (WC) profile, consisting of 74 
recording sites running along a WNW line (figure 1/2). 

ihe farst site was 400 km. from the shot. point and the 
last) 3600 km.” In addition, 9 sites made up an arc 1675 km. 
from the shot point. Stretching from Cold Lake, Alberta, 
to just south of the Trans Canada Highway near Suffield, 
Alberta, this arc encompassed approximately 20 degrees. 
Generally, the spacings between sites was 50 km., which 
was increased beyond the range of 2000 km. It was felt 
that this spacing would preserve the signal character from 
one site to the next, and this character could then be 
used to trace)ievents across the profile. Near Holden, 
Alberta, from 1754 km. to 1731 km., the spacing was 
reduced to the length of the array (2 miles). This 
continuous section of the profile was to provide data 
Overran anuicipavea Critical distance, {or varriveals from 
below a possible low velocity layer in the mantle, and was 
to be used for experiments in beam forming. In addition, 
two permanent sites were occupied for the duration of the 
experiment, one at the University of Alberta Edmonton 
Seismic Observatory near Leduc, Alberta (1840 km.), and 
the other at the intersection of the profile and the arc 
(WC 36). These stations were expected to supply data for 


stacking and for calibrating the spectral density estimates 
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Ofethe individual sources. 


The proposed recording sites were taken from 
large scale maps (4 miles = 1 inch) and located on isolated 
EW roads. It was the responsibility of the recording 
crews to make the final selection of the recording locations 
when they arrived on the proposed site on the day it was 
HOmbe cecupmed, in, choosing thesrine! site locations s-uhe 
crews had to consider such features as power lines and 
seismic background noise, while maintaining the proposed 
Spacings on the profile. Since these were small, it was 
considered best if each crew worked a particular range of 
the profile. In this way they could consider the 
final location of the previously occupied site in choosing 
the ew one, thtts-preserving "the planned spacings. The 
crews also had to estimate the recording times at each 
lecatiom, )Sinitially they sused tatvelocity -or 26. 2 km./sec. 
and an intercept time of +5 seconds to estimate the 
arrival times, as taken from a review of the Cumming and 
Kanasewich data. These times were to be modified and the 
recording gains adjusted by extrapolating the data collected 
during the exercise. On the other hand, crews occupying 
adjacent sites would have more information on which to 
base these decisions, and an experienced observer could 
review all the records and make these decisions while the 
crews slept. Also, one technician could then service all 


of the recording equipment. Ideally, of course, these two 
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Specialists would have no recording responsibilities. 


Both systems were employed. The three USGS crews 
occupied adjacent sites on the close range section of the 
profile (400 to 1000 km.), after which each crew took on 
profile assignments on the far section (2000 to 3600 km.) 
The University of Alberta crews were given profile assign- 
ments on the central section of the profile (1000 to 2000 
km, including the arc). Both systems worked well. It 
was found that the freedom of site selection offered by 
the profile section assignments was extremely valuable 
in the selection of good sites. However, due to the 
quality of data and the lack of time for adequate analysis, 
weak first arrivals were not always observed by the field 
parties, and hence no allowance was made for recording 
them on subsequent days. This type of oversight proved a 
handicap in the final data analysis, and could possibly 


have been avoided had the adjacent site method been used. 


IES) Field Equipment 


The USGS crews employed their standard mobile 
6 Station linear array refraction equipment. The Univer-— 
sity of Alberta crews used their mobile lz station refrac- 
tion y reflection linear array apparatus, the same as 
used by Maureau (1964) during his crustal refraction 
program and by Clowes (1966, 1969) during his deep crustal 


reflection programs. The latter arrays consited of le 
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Stations with spacings of 960 feet covering 2 miles. At 
station 6, as numbered from one end of the array, the leads 
of stations 2 and 4 were tapped, and a three component 
(vertical, radial, and transverse) station was established 
using Hall-Sears HS-10 seismometers. Elsewhere, other 
than at 2< and 4, Texas Instruments vertical S-36 seismo-— 
meters were employed, See ficures 1.3 and.) for the 
seismometer responses. All components of the array, in 
addition to timing information, were recorded on a RS-8U 
recording oscillograph. These records will hereafter be 
meterred tbo as the ficld monitor records. Ay Precision 
Instrument PS 207 AM/FM seven channel tape deck was used 

to record 6 components of the array in FM, and one channel 
of timing information in AM. Two modes of recording were 
available, one recording the odd numbered array components, 
hereafter referred to as the normal mode, and another, 

the reverse mode, recording the even components. Thus 
€ither a 6 station or a 4 stativon.(having in addition the 
two horizontal components) linear array, both with spacings 
of 1920 feet, were available for replay and digitization 


in the Edmonton laboratory. 


ee Field Operation 


The field operation, then, entailed the location 
of a noise-free site, setting up the array and the recor- 


ding equipment, and recording the shot. Noise records were 
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Pigure, 1.3 Frequency response of the Hall-Sears HS-10 


seismometers. 
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Figure 1.4 


Frequency response of the Texas Instruments 


vertical S-36 seismometers, 
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run to determine the state of the system and the background 
noise level at different filter frequencies. These recor- 
dings were run several times. Once before the system was 
completely set up, leaving ample time for repairs or site 
relocation, and again as soon as the station was in the 
ready state, in order to determine attenuations and filter 
Settings: for the recording of the shot. )Thesfirstshot 
was recorded in the normal mode, using the data collected 
at previously occupied sites as a guide for recording 
times and attenuations. The record quality and the 
eprival tamessof this) recordyweresthen checked-on the 
f2eld monitor record and used as a guide to vrecord the 
second shot of the day, which was recorded at the same 
Save.) lhe second shot. could, be» recorded Gn the reverse 
mode by reversing the array cable inputs; hence the name, 
reverse mode, Top priority was placed on the normal mode 
recording, as ateprovided the, largest array,on tape. If 
the first shot.record was of good quality,.the second 

shot. of the day, at: aiternate, sites. only, wasyrecorded 

in the reverse mode. A day's operation consisted of the 


recording of two shots on both tape and photographic 


paper. 


nee) Digitization 


The digitizing system consisted of a 6 channel 


hot wire recorder for monitoring the tape channels, a 
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two channel pen recorder for the tape timing channel and 
the digitizing time, aliasing filters, a 6 channel Packard 
Bedinimrloi plexer;varPackard ‘Bedus 1d) bat anadoe to Aaeatal 
converter, a clock (a square wave oscillator which governed 
the conversion rate), a digital counter, switching circuit 
between two Kennedy digital tape recorders, and the PS 207 
tape deck in the replay mode. Because the system was not 
designed to convert and record seismic data rapidly 

enough (over 100 samples per second (s/s) per channel 

werel requireds for statistical Significance, requiring 

600 s/s by the system), the FM tape which was recorded at 
30 inches per second (i/s) was replayed at 32 i/s. 

Thus the digitizing time was increased a factor of 8 

over real time, and a conversion rate of 90 s/s digitizing 
time was equivalent to 120 s/s per channel real time, 
which slightly exceeded the required rate. The aliasing 
Gael cer,s 24 dD, per octave, «corner frequency on 470 snZ 
digitizing time was used to assure that there was no 
signal of any significant power at a frequency higher than 
the aliasing frequency (7.5 HZ machine time). It should 
be noted that the data was recorded using a 24 db per 
octave filter with a corner frequency of 8 to 16 HZ which 
corresponds to 1 HZ to 2 HZ in digitizing time. Thus any 


signal near or above the aliasing frequency was severely 


attenuated. 


The data was digitized in channel sequence with 
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a flag marker added to the word of the sixth channel, and 
was recorded on the digital tape in this form. This was 
organized into records consisting of 10 blocks of 3600 
words each, giving a total of 36000 words per record or 
6000 words per channel (just over 50 seconds real time of 
seismic data). Imnter-record gaps are required between 
blocks. As the Kennedy recorders demanded more time to 
insert these than the sampling interval, two recorders 
and a switching circuit were used. When one recorder 
hadurecordéd a full) blocky ite put” invits!inter=record 

gap while the other recorder recorded the next block. 
This solution was considerably less expensive than using 
a high speed memory bank interfacing system. Thus each 


FM tape was digitized and recorded on two digital tapes. 


The data was taken in this form to the IBM 360 
Model 67 computer, and one tape containing all the data 
was built. During this process, longitudinal and horizontal 
parity checks and channel sequence checks, along with other 
IBM system checks, were made to insure that the digital 
data was reliable. On detection of an error the data was 
rejected and later redigitized. In this manner Lik digital 


records were compiled on two master tapes. 


1.6 Bookkeeping and Data Limitations 


Before data analysis could begin, a considerable 


amount of bookkeeping had to be done. The locations of the 
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shot points were taken from a technical letter published 
by the United States Geological Survey (Warren, Healy, 
Horiman, Kempe, Rauula, and Stuart; 1967), and the recording 
Sites were determined from large scale maps. From this, 
the profile distances were calculated using a Fortran 
program based on a Clarke spheroid model of the earth. 

The accuracy of the distances was limited to the scale of 
the map and was estimated to be 0.2 km. Time was synchro- 
nized between the shot point and the recording crews to an 
accuracy of 0.005 second using WWV as a standard. MIrre- 
gularities in the spreads were noted in the field log 
books, and the spread configurations calculated to the 
nearest 2 meters. The digital data was plotted on the 
CalComp model 663 plotter using the same time and amplitude 
scales as those of the field monitor records, Froma 
direct comparison of the two, the time of beginning and 

end of digitization was measured and the sampling interval 
calculated. Using this data, the timing of events could 


begin. 
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GHAPTER 23 


DATA ENHANCEMENT 


Prone IntEroCUuction 


Any signal consists of two parts, one being 
wanted and the other unwanted, the latter being defined 
as noise. The noise can be broken down into various 
Categories, depending on its characteristics; “in fact, the 
noise signal can be considered as a combination of several 
types. The superposition theorem tells us that we can 
treat the wanted signal and each type of noise as separate 
Signals. Signal enhancement, then, consists of the 
separation ofjithese components, or at least the suppression 


of the noise, 


Various processes have been derived for signal 
detection and noise rejection, These methods fall under 
two main headings, one being analog and the other digital. 
In general, for every analog process there exists an 
equivalent digital one, and vice versa. Thus this is not 
a true classification of procedures, but a classification 
of signal format and of the machine used to execute the 


procedure. 


It becomes a question of which machine to use. 


Analog computers require separate electronic components 
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for each operation, such as addition, multiplication, 
delays, etcetera, and wires for controlling the operation 
sequence, The machines tend to be expensive and very 
Specialized, with more noise being added to the signal by 
each amplifier. In digital computers, all operations are 
broken down into binary operations and the logic is 
represented by machine instructions generated by a language 
such. as Fortran,  Thus.one; machine: can.perform all) the 
operations, with the programmer varying the logic as 
required. Hence the digital computer is extremely versa- 
tile and free of noise. The only noise introduced is in 
the round off of numbers during the binary operations, 
which may be rendered insignificant. However, because of 
the flexibility of the computer and the break-down of 
Operations into simpler binary operations, the digital 
process can be substantially slower than the equivalent 
analog one. The selection of processes often breaks down 
imcoespeed versus flexibility, —For research purposes, 
flexibility is far more important, hence the analog signal 


was digitized and then processed on an IBM 360/67 computer. 


LAB Detection and Elimination of Spurious Transients 


Transients or spikes are characterized by their 
rapid change to a maximum followed by an over-damped return 
to, the normal signal level... These spikes occurred sduring 


field recording or during laboratory digitization. A 
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first difference method was used to detect the rapid 

change and a linear interpolation to replace the spike. 

A more sophisticated correction algorithm was not attempted 
due to the difficulty in establishing a suitable detection 
criterion. Moderate success was attained using the 
deviation of the first difference from the mean first 
difference as established over a two second interval in 

the record, this window being moved across the entire 
record. The allowable deviation was based on the standard 
deviation from the first difference. Once the spike is 
detected and its amplitude determined, and knowing the 

time constant and the damping of the system which generated 
the transient, it should be possible to calculate the 
transient signal. Thus it can be removed by subtracting 
the calculated signal from the seismic signal. It appears 
that it is possible to write an anti-spike algorithm, 


but the limitation of time prevented its full development. 


Ze) Frequency Filtering 


Local cultural activity and the wind generated 
background noise which was generally of a higher frequency 
than the wanted signal. A Butterworth zero phase band-pass 
filter, programmed by Alpaslan (1968), was used to reject 
this noise component. One filter with cut-off frequencies 
of 1.0 and 5.0 HZ passed nearly all of the signal and 


attenuated most of the noise. Another filter (0.5 and 3.0 
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HZ) slightly attenuated the Signal and rejected almost all 
of the noise; however, it also removed the character of the 
Signal. Figures 2.1 to 2.4 show the impulse and amplitude 
responses of these filters. Figure 2.5 shows some unfil- 
tered Early Rise records, while figures 2.6 and 2.7 show 


the large improvement after application of these filters. 


eal: Linear and Non-Linear Beamforming 


As this background noise was also incoherent,a 
method of beamforming which is based on this criterion for 
rejection is the simple delay and add routine. This is a 
linear method and amounts to timing the antenna (array) to 
Doinbeansa parcicular direction, Ihe acceptance lobe can 
be sharpened by using weighting coefficients on each of the 
Signals to be added. This procedure is commonly known as 
stacking or velocity filtering, as the delays depend on the 
geophone spacings and the assumed apparent velocity of the 
wanted arrival. Figure 2.8 shows the results of the before- 


mentioned records after application of this velocity filter. 


Another beamforming technique employed is a non- 
linear variation of the delay and sum. This method 
(Muirhead, 1968) consists of a delay, taking the N” root 
of the signal, summing, and then raising the result to the 
Nn” power. Hereafter it will be referred to as the N* root’ 


stack. As can be seen in figure 2.9, this technique had 
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Figure 2. 


Impulse response in the time domain for a 
four pole zero phase Butterworth 1.0 - 5.0 


HZ band pass filter. 
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Impulse response in the frequency domain 
for a four pole zero phase Butterworth 


0 = 5.0 HZ bandpass filter. 
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Pigune 2.3 


Impulse response in the time domain for a 
four pole zero phase Butterworth 0.5 - 3.0 
HZ band pass filter. 
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Figure 2.6 


These same records after applying the 


U0 War 5,0 HZs Butterwontie tGer, 
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Figure 2.8 


These same records after applying the 


Velocity filter: 
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Figure 2.9 These same records after applying the 


a root stack filter. 
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ereaver NO1se rejecting ability than the normal suack, and 
also haa better velocity resolution, It achieves Uhis at 
wie expense Of Signal distortion, Une, rejeculonand 


GaSvOrulOn Increasing as a puncuLon,o1 IN. 


The outputs of the various stack routines were 
displayed using the CalComp’ 663, digital “plotter® © Here each 
output trace is the result of assuming different spread 
velocities. The sensitivity of the beamforming filter is 
limited by the long wavelegths (high velocity, low 
frequency) and the shortness of the spread. A geophone 
separation of one or two kilometers would have been more 
appropriate, as seen in retrospect. However, the results 


were satisfactory, as verified by the figures. 
Let/us consider the stacking procedures in more 
detail. Consider K channels of data, 


X, (t) jet, K en 


with channel j having gain G. Let the channels be sampled 
at m samples per second for T seconds, giving us a K by M 


maori 


where M is the total number of samples per channel and is 


given by 
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For the N” root stack, the intermediate result is 
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Note how the sign is conserved in the N root (equation 


2.8) and the N” power (equation 2.9) operations. 


Let us consider the suppression of a spike of 
amplitude A at i=I on channel J. Assuming all signals 
are the same except for phase, we have, letting Cj=L 
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represented by a Gaussian distribution of spikes, it can 
be shown that the signal to noise ratio is improved by 


apeactor of #7 (Smith, 1950). 
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Therefore spike rejection varies as ane. However, in using 
the crude approximation (equation 2.17), the filter 
operator was changed into a linear one, whereas, in fact, 
the filter is highly non-linear. A much more sophisticated 
treatment is necessary to gain any more than the general 
behaviour of these filters as presented here and in 


figure: 2.10. 


Za Display of Results 


The display is just as important as the elimina- 
tion of noise, because the organization of the display 


plays a major role in signal interpretation. All displays 
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Figure 2,10 


Spike attenuation using the NY root stack 


on K channels of input. 
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were generated on the IBM 360/67 computer and plotted 
using a CalComp model 663 digital plotter. Even the 
unfiltered records could be improved by organizing the 
record display, so that the center axis of each trace was 
displaced from the time axis a distance proportional to 
the displacement of the corresponding geophone from the 
center of the spread; hence the apparent spread velocity 
could be réad directly from the display. This proportional 
factor and the gain level could be adjusted to ensure that 
each trace was distinct from another and event character, 
frequency, and amplitude could readily be compared on all 
records, This was not necessarily the case with the field 
monitor records, particularly when the recording was made 
in the reverse mode; in this case the traces were 
scrambled, and it required the flexibility of the CalComp 


plotter to display these records in a meaningful manner. 


Signal interpretation involves the identification 
and timing of events, which could not be done on the 
computer. An experienced observer, through a visual 
display, is able to perform complex and poorly understood 
(therefore unprogrammable ) correlation functions to 
identify these events and trace them across the profile, 
provided that the signal to noise ratio is large enough. 
The enhancement techniques are used to unmask the signal 


from the noise, but cannot by themselves interpret the 
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data; hence the importance of the visual displays. They 


are the tools of the interpreter. 


Figures 2.5 to,2.9 aré characteristic of the 
displays Of individual records, To give a picture of the 
entire profile, a reduced travel time composite plot 
(figures 2.11 and 2.12) was composed using the best 
recordings, which were band-pass filtered and stacked. 
These and similar displays of the other records were then 


used in analysing the data. 
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Haeurers2.i1 Composite plot of velocity filtered records. 
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CHAPTER 3. 
DATA INTERPRETATION 


ok Earth Model - Ray Theory 


Gonsider a radially symmetric earth made up 
of infinitesimally thin concentric homogeneous shells 
through which we shall pass a ray obeying Snell's law. 
Denote the physical properties of the i? veheld (between 
Tadka ies ard r;) with the subscript i; hence the angle 
of incidence and the angle of refraction at the to 
interface can be denoted as a; and at,, respectively, 
the angle being defined as the angle between the ray and 
the normal to the interface. Further, let the intersection 
of the interface and the ray be denoted as 5S; (see figure 
Su lleeetegiven ray, then, can bevderined by the angle of 


incidence at the source, a,, at radius r, or depth r.-r.. 
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Breure 3.1 


Diagram of a ray passing through an earth 


made up of homogeneous concentric shells. 
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where QO is the perpendicular distance from the ray segment 


in the i* shell to the center of the earth (figure 3.1). 


Therefore, Snell's law becomes 
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where p is a constant known as the ray parameter, which 
is unique for each ray. Now we can drop the subscripts 
and consider the velocity function v as a continuous 
muneeLon Of ry vhav 2S, Vv = v(r). Then the angle (ee 
incidence is also a continuous function of r, (an rane) ys 
and the bottom of the ray, or the minimum radius OL) Ulie 


ray, ry, is defined by 
a.(r,) = a /2 BO 


Hence the depth of deepest penetration may be defined 


as Yo-r,- Then 
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That is, the ray parameter is defined by the radius of 
the source and the angle of incidence at the source, 
and is a constant everywhere on the ray. Then by the 


bottom condition stated in equation 3.6 
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Equation 3.10 is in the form of a transcendental equation 
which, depending on the function v(r), may or may not 
be solved explicitly for rn. In cases where it cannot 


be solved explicitly, it may be solved by iteration. 


Consider the ray path S originating at the 
polar co-ordinates ge 5x0) and with an initial angle of 
incidence a,. An element ds of the ray path can be 


expressed as 


ds = dre 4% + Hem is ok 


where u, and Ug are unit vectors in the r,@ directions 
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respectvely. Two identities are immediately apparent 


and are needed in the following development. 
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(ds) = (dr) + r (Ae) 
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from the definition of p given in equation 3.7. 


substituting equation 3.15 into equation 3.12 
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Therefore the angular component of the ray path from its 


source to its depth of deepest penetration 1s 2. ven! by 
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Solving for'de an equation 3.1359 substi Guvineg 
iG into equation 3.12, and using! equatiom 3\.//-v0 eliminate 


a(r), we get 
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The time dt taken to traverse the path lengthyds is 
given by 
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Therefore the time taken to traverse the ray path from 


the source to the depth of.deepest penetration is given 
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For a complete ray path, from source to receiver, 
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where r. is the radius 6f the receiver. Whenwc. = 2. 
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this simplifies to the form 
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The equations of @ and T in the forms of equations 
pyeoeend jecOr are the most useiul, isanee dattercny 
velocity functions may be used at different angular 
displacements from the source. That is, v(r) can be 
varied as a function of @ and becomes v(r,@). However, 
unless the ray parameter p is also varied, it LSS Shall all 
assumed that the shells are radially symmetrical, but 
that the radius of the shell can be made to vary With oO 
for individual ray paths. In other words, tite eid ot 
the shell has not been taken into account unless the ray 
parameter p is suitably varied, Hence the usefulness of 
equations Bae Jane 3.30 51s limited to varying the elevation 
and crustal structure under the source ana the receiver, 
still assuming beds parallel to the surface of the earth. 
This can be accomplished by using the source velocity 
function in the down integral (tr. Go r,) and using a 
suitable velocity function in the up integral as based 
on the expected angular displacement from the source 

to the receiver end of the ray, as estimated by 2Q, 
from equation 3.3¢. It should be possible to consider 


an asymmetric earth by varying the ray parameter p as 
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a function of the dip of the velocity shells as represented 
by sherical shells displaced from the center of the 
earth, 
ae Solutions of the Ray Path Integrals 


Note that both @ and T are of the form 
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There are several ways of evaluating this integral. 


The obvious method is to evaluate 
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However, f(x), or rather f,(r) and f(r), are not well 
known functions, as v(r) is assumed at discrete values 

of r, and some form of extrapolation is used between 
these points. Therefore vir) is, at best, a continuous 
function with discontinuous high order derivatives. Also, 


integral solutions exist for only a limited number of 


\ 

7 a ne wy a . | 
att aon ‘eifente’ xo 30 av (3a ie ‘a % 
re 00 “toToeo aie ie 7 oi es q ne ai ts Le 


Sot 
[= 
! 


a 


a? ity hy i 
0) ae 
alstsedal ated yeh oat To dso tisha” Bok 


miot eft to eis T bse YS atod weds oto 


~Leraotai eindt aniteulava to sysw istevee ote ore 2 


steavievea ot sf bodega evo byaia ae T 
ie c Ma lo 4 | ub » = 
. - » (+) * = \* } 1 


ef 6€ (A) A) =) Ay, in 


fiew ton ots : (1). bs: (t)g? “xsittat xo t (RYE. of ¢ = c 


~ 


zeulsv stetoeib Ge boywiees ae tt} 258 , ato b tome sn mm T: 


ares boath at hobaski sqentes Yo mo (omer ‘he ere 


oh avo ones 8 ods: até Yeu (k)y swine ei ee 
: 93th Pans ‘ealoxin ae ons ne i trod 


. 2° home 7 er p Q 
. 3 _ - dak 5 - re oe ee a ere » ba T 
ne Flere is ia ae ee ot Tap Ju! ; 


\ vs) 


7 : a 7 
7 ' 
i : 


L7 


velocity functions, and these may not necessarily be 


poivsmcaliv realistic, 


Let use-consider a. finite number of concentric 
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travelling tangential to the curve of r-=-r,. Ditsacan 


be considered as a special case of reflection where the 


Secle Gieancidence is a/2. 


Let us consider the solution to the travel 
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cannot be readily integrated, and no solution has been 


derived here. 


As expressed in equations 3.18 and 3.26, both 
Ueand Tare of the form. Ko= ps f(x)dx. When the function 
v(r) is not well known and various methods of extrapolation 
between radii of assumed velocities are used, then a 


center point integration formula of the form 
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Using this method, the integration of £(x) Ws equated to 
thesarea Under the curve f(x), and this is approximated 
by the sum of the areas of several rectangles of height 
Ge )> 2 = 0,N “and of width ax) Both positive and 
negative errors occur which tend to cancel each other 
(see figure 3.2), and which are further reduced by 
decreasing 4x. The accuracy is then limited by the 
accumulated round off errors, which increase as N is 


increased and hence as Ax is decreased. 


With a knowledge of the behaviour of the 
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Way be caken with small sacritice to Uhe accuracy sor ‘the 
approximation and with a large decrease in the accumulated 
round off error. As the pole is approached, the curvature 
Sf obhe Luncb1on increases, Hence, smaller ssteps arc used, 
and the error of the approximation is reduced, Also, the 
round off error is small, because each of the terms 
representing the area under the curve is of small width 
but of large height, therefore of considerable area. 
Another range of x exists, where the function is of 
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Wagure 3.2 Center point integration formula approximation 
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when the pole exists just below a low velocity shell, in 
which case there is a region immediately above this shell 
where finer increments in the integration are required. 
Ghuse by Suitably varying Wx alt ds spossi blero mi namize 


both the error of the approximation and of the round off. 


Fortran programs were written (see appendix) 
using both the thin shell constant velocity and the 
center point integration formaulae. To check the accuracy 
of the programs, they were used to reproduce the travel 
times of Jeffreys and Bullen(1958); Lewis and Meyer 
(1967); Johnson (1967); and Herrin, Tucker, Taggart, 
Gordon, and Lobdell (1968). Furthermore, the two 
methods closely agreed with one another in each of the 
above checks. In addition, the center point integration 
formula agreed closely with the exact solution Ol “une 
constant velocity formula in the case of a homogeneous 
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Se: Ray Flux Calculations 


In conjunction with the travel time curves, 
ray paths and amplitude predictions were made. The ray 
paths were calculated by evaluating Bo toateeach step) or 
the integral and the results were plotted using the 
CalComp digital plotter. The energy density was calculated 


from the ray flux, which was shown by Bullen (1963) to be 
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where I is a constant of proportionality denoting the 
energy per solid angle emitted by the source, a, is the 
angle of incidence at the source, a, is the angle of 
emergence at the receiver, r, is the radius of the 
receiver, @ is the angular distance from source to 
receiver, and da,, d@& are increments imea,,6 respectively, 
To derive this equation we consider all the rays leaving 
the source between the angles of incidence a, and a,tda, 
(see figure 3.3). The energy radiated in this ray 

bundles 2s 2cri-sin a, | das], the modulus of da, being 

used to insure that this term is positive. This ray 
bundle arrives at the receiver radius r, between the 
distances 6 and @+d@, which sweeps out an area On a 
sphere of radius r, given by 2nr. sine|de]. Again, 

the modulus of d@ is used, as d@ may be negative, whereas 
the area is always positive. Then the cross-sectional 
area of the ray bundle is given by Qtr, sine Cos a,|do| 


and the energy density at the angular distance O7is 
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Fagure 3.3 


Diagram of two rays passing through the 
earth, showing the change in the angular 
displacement of the emerging ray Lor a 


small change in the angle of incidence. 
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Finally, the square of the amplitude is proportional to 


the energy density, whereby 
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Using angles of incidence a; , a,tda., and a,-da,, rays 
were passed, which arrived at the receiver radius with 
am angular distance of G, &., and G.. respectively. 
Hence d@ could be estimated using the first difference 
method on either side of a, and smoothed by taking the 
average of these terms. Thus all the arguments of 
equation 3.77 are known and A(@), along with the travel 
time curve, can be calculated for any given velocity 
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Bie Interpretation of the Travel Time Curve 


The various earth models were derived by 
varying the velocity function until a good fit was made 
between the theoretical travel time curves and the 
events picked from the records. Model U of A #1 was fit 


to the first arrivals and to a prominent secondary 
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arrival as reported in the Technical Letter #6 (Warren 
et al., 1967). The ray and travel time diagrams are 
shown in figures 3.4 and 3.5, the latter also showing 
the Lewis and Meyer (1967) model for comparison. The 
Lewis and Meyer profile was also a part of Project 
Karly Rise and was run along the (Weis parallel (see 
figure 1.1), which is close to the Western Canada 
Doone.) [heret ere, the two prot ples .siould snow 
similar results, 9as. im fact, theysdO, as Ullustrated: im 


the diagram. 


When a full analysis of the secondary arrivals 
had been made, it became evident that a simple model 
such as a modification of U of A #1 could not explain 
all-of the arrivats: -Invordér to ‘explain the shadow 
Zone and Cuspingein thestravel Urme “curve, one needs 
velocity reversals and large velocity gradients. Model 
U ot Ao3 produced the best fit to the observed data 
(figures 3.4 and 3.6). This model is markedly differnt 
from U of A #1. By removing the sharp corners of the 
velocity function of U of A #3, a modification, U of A #3a, 
resulted. As seen in figure 3.6, there was no significant 
difference in the resulting travel time curves of these 
two models. When the ray criterion of explaining events 
was relaxed, and the possibility of diffraction was used 
to explain the extreme ends of the cusps, then U of A #2 


produced the best fit to all parts of the travel time curve 
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Velocity and travel tame curves sor vue 
models U of A #2, U of A #3, and U of A #3a 
along with the Western Canada profile 


obSservaclons. 
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excepy the vips Of The Cusvs, as Shown by Dieures oo end 
3.7. The selection of a unique model to fit the data and 
possibly the determination of the shape of the transition 
zones lies in the matching the predicted amplitudes to the 
observed ones, as well as maintaining the fit to the 
observed travel time curve. This work is presently in 


progress. 


Figure 3.8 is presented to summarize this 
work. Lo shows. U of Ao, Which Gs a imodel t1tcing 4 
continental profile, and compares it to a model worked 
out by Johnson(1967), which is based on an oceanic 
profile. Also, a model based on a summary of the United 
States continental data (Herrin et al., 1968) is shown 
for further comparison. The rough features of our model 
and that of Johnson are seen to coincide, with ours 
showing a delay of about 2 seconds. The Herrin model 
lies siiently above modell U ef A #3 and has none of the 
detailed cusping of our model. Further work in this 
area should clarify the details of the upper mantle, 
but the gross features must be as represented here. In 
panticular, a reversal of this prefite should clarify 


the true velocities of the various zones. 
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Figure 3.7/ 


Ray diagram for body waves passing through 


the earth model U of A #2. 
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Veloerty and travel time curves) on the 
final model produced by the University of 
Alberta (U of A #3) along with those 
resulting from the oceanic model produced 
by Johnson and the average crustal model 
produced by Herrin, Tucker, Taggart, Gordon, 


and Loodell. 
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APPENDIX, 


PORTRAN LV PROGRAM FOR THE TIME AND DISTANCE INTEGRALS 


The following Fortran IV program was used to 
calculate the travel time curve for a given earth model. 
It integrates from R, (which isjequivalent to r, in this 
Hepen) toa, calculated revusing velocity model y(rj0)},5 and 
thenshas the option of resetting the velocity structure for 
the up integral by evaluating the v(r,@), where @ is 


twice 6. ‘ 


Two*integral formulae were used; ana hence the 
listing of two down integral and two up integral subrou- 
tines. One is the center point integration formula and 
Ghe Other is “the thin homogeneous Shell. AGhird: up 
invegral’ subrowtine, which simply doubles"the value” or 
Sere and T.,, is also listed. This is valid only when 
v = v(r) and r,=r,. The appropriate subroutines must be 


selected by the user. 
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Ga Gay (Gauc) YEG) AG Gat ee? ry any SD GON aR 


Signa e 


N.2 


REAL*4 RDEL(50,20),VDEL (50,20) yDELMOD(50),RMAN(B1),VMAN( 81) 

REACES TR ULOL) sVGLOD) 

RFAL*8 RO(101),VO(101L) 

RFAL*®8 RADDEG RE R1LyR2yRZ eVREDeTINC AT LyAI2,AINC DAL yAL yDEL TA, TIME 
REAL*®B NEL TAC yDELTAMyDELT AP, EM, EPy SLOPE yP RAD, DELTAKy DELTADyTRED 
REAL *8 V1lyPE2,TESTyDRyRADC »RADTy VELCe VEL TyRDVeAyTHOLD ey VEL oT IMEC 

100 FORMAT(5F10.0,3110) 

103 FORMAT(*OANGLE OF INCIDENCE = *,F10.5) 

LOGOFORMAT("1LRADIUS OF THE EARTH ="-Fl10.39'KM*/*"OSTART INTEGRATION AT® 
LePLOs3e2KM,tSTOP ATESFLOS3SSUKMEMRADIUS OF PENETRATION LIMITED 70%, 
2EL0.3,"KMY/Z* L/DR STARTS AT'yI4y" STOPS AT',18,* & ITEST =*,18) 

1050 FORMAT(*ORAY FLUX GN EACH SIDE OF THE RAY TS %42(D12.6% 

Ly = UGG Pe yeelors sti 85" ~AV FUUK =", F10e5) 

106 FORMAT(*1 "3 

107 FORMAT(*LRAY PATH COORDINATES DELETED") 

108 FORMAT(2F10.0515,54X,11) 

LLOOFORMAT (*=",20X,"RAY STARTS" yFBe29" DEG "95Xy*RAY ENDS* FB oC, 

LA BWEGS Es SKF? INTERVAL YNOS 4153 
J$=50 

1$=20 

K$=81 

L$=101 
RANDEG=3.1415926536D0/180.0D0 


DATA TGARD ML ELOSOCNRAY 
NRAY IS THE NUMBER OF RAYS PASSED PER SHELL FOR TRAVEL TIME CURVE 
NRAY T= 1 IF. UNSPECIFIED (BLANK FIELD) 


READ(5,100)XOUT 
NRAY=IFIX(XOUT) 
TF(NRAY.LT.2)NRAY=2 


BATA TCARDSZ24ELOsOL XOUT 
IF XOUT NOT SPECIFIED, THE RAY PATH CO-ORDINATES ARE DEVERED 
TF XQUT NE 0.0 THE RAY PATH CO-ORDINATES ARE CALCULATED 


READ(5,100)xX0UT 
TFC LTEIX(XOUT) -EQ.0WRITE( 6,107) 


DATA CARD) 3S,FLOSO°RE 
REIS THESRADIUS OFETHE SPHERE 
RE = 6371-0 (KM) IF UNSPECIFIED 


READ(5,1LOOIRE 
LF (RE 2EQ.0.000)RE=6371.000 
2 CONTINUE 


DATA CARDS SET 1 TITLE AND DATA OF VELOCITY STRUCTURE 
SEE SUBROUTINE DATA 


CALL DATA(RDEL »VDEL yDELMODs J$¥1 $e RMANyVMANy KS yRoV LS) 
1 CONTINUE 


OATARCARD SET: L-+ 4 5F10.0,3110 R1yeR2eR3BeVREDeTINC,IDRL» IOR2,1 TEST 
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Ripe or RACLTUS (KM) GF SOURCE 

RE TSORADTUS? (KM)? GEM SPHERE 

Ris eater tiryUNSPECTRIED 

PBrIS RADLUSOAKM? OF RECEIVER 

BeeeakeEDTEAUNSPEC TELE 

R2 TS MINIMUM RADIUS (KM) RAY MAY PENETRATE 

VRED IS REDUCTION VELQCITY (KM/SEC) FOR REDUCED TRAVEL TIME 
VRED= SSS. 5SKM/SEC ULF UNSPECEFIED 

ORAM TSOTORRAIGSTART OF INTEGRAL 

WHERE DR = L/TOR , IS THE INCREMENT IN THE INTEGRAL 

[DRIP ey2. 1, UNSPECIFIED 

[LOR2ZVOS*IOR=LIMIT FOR END OF INTEGRAL 

LOR27=f2E28N1F UNSPECIFIED 

ITEST IS NUMBER OF INCREMEATS TO BCTTOM OF RAY (REND) BEFORE 
IDR = IDR*2 

BRESTO=) LOOPTECUNSPECLELED 


READ(5,1LOOIR1»R2,R3y,VREDyTINC, IORI IOR2,1TEST 
LE CRISEO. OLODOIRL=RE 

TF(R3.EQ.020D0)R3=RE 

IF (VREDJEQ.0.000) VRED=8.5D0 
LFCTINCSEQL 0.000) TINC=0.5D00 

TFC IDR1.E950) IDRL=2 

IF( 1DR2.EQ.0)10R2=128 

PEGIPTESTSEGSOVETESTHLO0 
WRITE(6,104IRE,R1gR3eR2yIDR1Le IOR2, TEST 
CONTINUE 


DATA CARD SET 1 # 5 2F1000515954X,1T1 ALL,AT2sNRAYS,ISTOP 
PRIS TOP SUP IERMINATE PROGRAM 

Leeis.OP 8 , RE-ENTER PROGRAM AT DATA CARD SET 1 + 4 

IF ISTOP = 7 » RE-ENTER PROGRAM AT DATA CARD SET 1 

AIL 1S MAXIMUM ANGLE OF INCIDENCE (DEGREES) 

AL2 IS MINIMUM ANGLE OF INCIDENCE (OEGREES) 

IF AIL IS UNSPECIFIED, AI1 = 89.59 Af2 = 20-0 

NRAYS IS NUMBER OF RAYS PASSED PER SHELL FOR FLUX CALCULATION 
NRAYS = NRAY IF NRAYS.LE.NRAY AND RAY FLUX IS DELETED 
KEKEKRAY FLUX CALCULATED CNLY IF NRAYS.GT .NRAY €# EERE EH 


ot 


READ(5,lLOB)AITL,AI29NRAYSoISTOP 
TELTSTOP.ERSOICALL EXIT 
IF(ISTOP.EQ.8)GO TO 1 

Ret Toque EO eT GO 10:2 

IF (AI1.GT~90.0D0)AI1=90.0D0 
{F(AT1)1010,1010,1011 

CONTINUE 

APT 7895 500 

AI 2=20.0D0 

CONTINUE 

AT1=AI1L*RADDEG 

AIT2=AlL2*RADDEG 

LF (NRAYSoL TeNRAY )NRAYS=NRAY 
GALL SETIR CEL pVDEL pb $9 1 So RMANy VMANGK$ pROeVOeL Sole IDIMO) 
EDUMOM S=CL0IMC P= 1 
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LOl<2 


INTR = 1 
PU LLOOGe Tes 25° LO LMOM 
AT=GROCE) ¥VOCLIIS(VOCT)I*RO(1)) 
IF(AT.GT.21-0D0IGN TO 1000 
AT=DARSIN(AT) 
LeCcaAr.GT.ALRLYGO EG "1000 
LEUCAL LU WoARAVAT=AALZ 
AINC=(ATI—ALI)/OFLOAT(NRAYS=~1) 
DAT=0.1L0DN0*RADDEG 
IF (AINC.LT DAT JOATHAINC 
IF (NRAYSJEQ.sNRAYICAIT=0.0D0 
LARA T=} 
PA PCINT SRS SNES INTRVL) GO TO LOOT 
ANGLEL=ATL/RANDEG 
ANGLE 2=AI/RADDEG 
WRITE(6yL1LOYANGLEL sANGLE2,INTRVL 
ODO 1001 IT=lLyNRAYS 
TA=!II-1 
AIL=AIT1-AINC*¥DFLOAT(IA) 
AL=AI /RADDEG 
WRITE(6,103)Al1 
AT=AILL-AINC*¥DFLOAT(IA) 
Eel 1 .eGel JAI=Al—1.00—8 
[IFC ITI-EQeNRAYS JAT=AI4+1.00~8 
IQUT=LFIX(XOUT) 
OCALL RAY (RDEL »VDEL sDELMOD yJS_1$¢RMANYVMANG KS pKa Vol SeAT RI gR29R3¢ 
LVRED, TINCT OUTs IDR1, IOR2, 1 TEST» RADDEG, DELTA yT IMEyREyP) 
NELTAC=DELTA 
Tee C= eM E 
TE(DAT.E0.0.000)G0 10 5 
1@=1 
OCALL D VEL{(RyVylSs 
123,V1l,IC,10IM) 


VL=ANGLE OF EMERGENCE 


VL=DARSIN(P¥V1/R3) 

DELTAC=DELTA 

AT=AL1-AINCXDFLOAT(ITA)-DAI 

1OUT=0 
OCALL RAY (RDELeVDELyDELMOD J$y1 $¢RMANG VMANS KS yRo Vol Se AT yR1 yR29R3y 
LVREDsTINCy LOUT, IDRLe IDR2y ITESTs RADDEG »DELTAyT IMEgREyP) 
DELTAM=DELTA 

Pel ietOs1)60 70° 10)2 

AT=AI1-AINC#DFLOAT(IA)4DAI 
OCALL RAY (RDELsVDELyDELMODs J$p1$+RMANyVMANGKS pRoVoLSeAT Rly R29R3y 
LVPED»e TING, LOUTs IDR1e LDR2y ITESTs RADDEG »DELTAsT IMEVRE oP ) 

CONTINUE 

DELTAP=DELTA 

EP=DABS (DAL¥RADDEG) DS INCA TL *RADDEG) /0COS(V1) 

EM=EP/LOSIN(( DELTACH0ELTAMI/ 2.000) *DABS ( DELTAC-DELTAN)) 

EP=EP/(DSIN(( DELTAC4DELTAP )/2-000) *DABS( DEL TAC-DELTAP)) 

EMLOG=EM 

EMLOG=ALOG1O(EMLOG) 
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1000 


EPLOG=EP 
FPLOUG=ALOGLO(EPLIG) 
AIL=AL1L-AINC¥DFLOAT(TIA) 
DIST=DFLTAC#RE 
TRED=UTIMEC=DELTAC*RE/VRED) 
LeU eS LIVEPLOG=EMLOEG 
ITF(TI~EQs.NRAYS) EMLCG=EPLOG 
FLUX=(EPLOG+EMLOGI/2.0 
WRITE(6,105)0EP,EPLOGyEMyEMLOGs FLUX 
MODEL1L=1 

CONTINUE 

CONTINUE 

INTR = INTRVL #1 

AT1=Al 

CONTINUE 

GO TO 3 

END 
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OSUBFUUTINE RAY(RDEL,VLCEL,DELMOD,J$y 1 S$yRMANs VMANoK Seo Vol Sy 
TAL sR1yR2,R3,gVRED,TINC, LOUTs IOR1,IDR2, ITESTsRADDEG DELTA, TIME, REsP) 


CALCULATES TIME, DISTANCE AND RAY PATH FOR A RAY OF GIVEN ANGLE 
DE eINGCEDOENCE 


REAL*4 RDEL (56516) sVOEL(S$,1$) 2,DELMOD(S$)9RMAN(K$) > VMAN(K S$) 
REAL*8 R(L$),VIL$) 
REAL *& RPAODEG Rey Ry R22 R3,VREDr TING ALI; Al2,AING DAT ATS DELTA, TIME 
REAL*8 DEL TAC »DELTAM,DELTAPy EMyEPy SLOPE yPyRAD,DELT AKy DELTAD,TRED 
REAL*& V1l,PFE2,TESTyDRyRADC sRADT,VELCy VELT pROVeAsTHOLD, VEL, REND 
LOOGCFORMAT(*O® »5X,* DELTA =*-F8.2_"KM T-=DELTA/',F4.29* = "yFle2y 
Eee SECY TRAD VOIMr =", F lacs KM?) 
LOZOFORMAT (8 4+" 478X—9 "DELTA ="yFBe2e"KM T=DELTA/'yF4e20" = "y 
Dewy eee o wo 
MODELIL=1 
OLF (MODEL .NEsMODELIICALL SET(RDEL »VDEL yJ$_ 1 S$pRMANy VMANGK SyRo Vol $y 
YMODEL1,IDIM) 
MODEL=1 
P=0.000 
DELTA=0.0D0 
TIME=0.000 
THOLD=0.0D0 
IFCAI.GE.90-0D00*®RADDEG) RETURN 
OCALL NOWN IT(RyVol So 
LAL eR1l_yR29R3 eI DR1y IDRZ2,1TESTe I DIM, LOUT,DELTAsTIMEyREePyRAD, THOLD, 
2TINC,IDReIPEFLT,REND) 
DELTAK=RE®DELTA 
NELTAN=DELTA/RADDEG 
TRED=TIME-DELTAK/VRED 
WRITE (6,100) DELTAKsVRED»eTREC,RAL 
DELTAK=DELTAK+4DELTAK 
Ic=1 
OCALL D VEL (RV yl $y 
LREND,PE2,!1C,IDIM) 
IF (IREFLTeNE.~O)PE2=DARSIN(P*PE2/REND) 
OCAtL MOD(DELMOD,J $y 
1DELTAK »MODEL) 
OILF(MODELeNE-MODELLICALL SET(RDEL» VDEL »J$e1$yRMANyVMANo KSeR a Vol So 
1MODEL,IDIM) 
2 CONTINUE 
IE (MODEL ~EQeMNDEL1)GO TO 5 
PG=1] 
OGALE (D VEEURy Vat o> 
LRENDsVleIC,IDIM) 
IF CIREFLT.EQ.0)P=REND/VI 
IFC IREFLT~NE2O)P=(REND/V1) ¥DSIN(PE2) 
5 CONTINUE 
ITOR=1 
DO 3 1l=l,1DR2 
INDR=IDR*¥2 
Pe((DRsGeslDR2IGO 10:4 
3 CONTINUE 
4 CONTINUE 
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1 UCONTPINUE 
ORES TPE Ei CRs V 51S 
LAT gR1yR2,R3,LORL, IDR2,1 TEST» IDIM, IOUT» DELTA,TIME,REsPyRADy THOL Dy 
2PTINC s FOR, IREFLT sREND) 
TDR=1OR1 
RAN=R 1] 
DELTAK=RE*DELTA 
DEL TAD=DELTA/RADDEG 
TRED=TIME-DELTAK/VR&ED 
WRITE(65102) ODELTAK,VRED, TRED 
RETURN 
END 
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OSUBROUTINE DOWN I(R Vel $> 
TAL» RVGRZ2.R3s1DR1, ILORZ, 1 TEST» IDIM, IOUT,DELTA,TIMEsPEeP »RADy THOLD, 
21 INC sy hOR,IREFLT»RENDI 


INTEGRATES FROM RL TN REND CRAY BOTTOM) 
USING CENTER POINT INTEGRATION FORMULA 


REALS R(LS)sVILS) 
QREAL*8 RADDEG sREgR1 R29 R3 eVREDe TINCe ATL eALZ,AINC DAI yA,» DELTAs TIME 
REAL *8 DEL TACsDELTAMyDELTAP,EMyEPy SLOPE 9 P oR AD, DEL TAKs DELTADsT RED 
REAL *8 V1 yPE2e TEST» DRyRADC yRADTs VELC p VELT eROVe Ay THOLDy VEL» REND 
100 FORMAT (LX s50 7X 9X8 7Xet* Zt yg TXe* Th) 
101 FORMAT(8(1X%,014.8)) 
Ic=1 
TN=0 
ITF (P1L.EQ.0.0D0)IR1=R(1) 
OCALL D VEL(RyVelS, 
IR1,Vly,ICyIDIM) 
IF (P2EQ.0.0D0)P=H=(RI/VILI*DSIN(CAT) 
PE2Z=P*P 
CALL BOTTOM(R,V,IDIM,P,AI,REND,IREFLT) 
RAD=R1 
IT=IC 
TESH=DEPLOA TAT TEST? 
TOR=I10R1 
DR=1/7DELGAT(IDR) 
IF (IOUT. EQ.) WRITE( 6,100) 
4 CONTINUE 
RADC=RAD-DR/?2 .000 
RADT=RAD-—DR*TEST 
OGCAUE FD VELCRSVsCS, 
LRAODC,VELC,IC,IDIM) 
OCALL 1D VEL(R»Vel$e 
LRADT,VELT,IT, IDI) 
T[EURADT/VELT.-LEsPIGU TC 1 
TEURAUT/ VELOC. EE sPIG0 10 2 
IFCRAD-OR.LT.R2IGC TO 3 
ROV=RADC/VELC 
A=RNV¥RDV-PE2 
A=DR/DSQRT(RDV*RDV=PE2) 
NELTA=DELT A+P *A/RACC 
TIME=TIME+RDV¥A/VELC 
RAD=RAD-OR 
IF (TOUTLEQ-1) CALL OUT (DELTA, TIME,R Ee TINCy THOLD yRAD GID 50) 
GO) TO 4 
i GON TINUE 
2 CONTINUE 
IT=1 
LDR=1DR*2 
DR=1/D0FLOAT(IOR) 
PEGTEST~£001-000)GD 1G 5 
IF(IOR~GE~IDOR2)TEST=1-0090 
GO TO 4 
5 CONTINUE 
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[F{ TOUTS EQe1LICALL CUT( DELTA, TIME eREsTINCs THOLDy,RADyIDy1) 
RETURN 
CONTINUE 


MODIFY LAST STEP IN INTEGRATION TO THIN HOMOGENEOUS SHELL FORMULA 
USER THIS FORMULA OMY TF  TREFLI =0 
OTHERWISE RETAIN PRESENT FORM 


DR=RAD=R2 
RADC=RAD-OR/2 2.009 
OCALL OVEL(R Vel, 
LRADC,yVELC,IC,IDIM) 
RDV=RACC/VELC 
IF (RDV*ROV-PE2.LE.0-0D0IGO TO 6 
A=DR/DSORT(RDV*RDV=PE2 ) 
DELTA=DELT A+P*A/RATCC 
TIME=TIME*#RDV¥A/VELC 
RAND=RAD=—ODR 
CONTINUE 
[FUL OUTCEQS TIC ALL QUT (DELTA, TIMEsREyTINCs THOLDsRAD,IDs1) 
RETURN 
END) 
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OSUBROUTINE UP YT(R VylL hye 
LAT RY «R2;83,1CR1, FORZ, LTEST s+ 101M, 1DUT, OELTAsT IME, REsP y RAD, THOLD, 
2T INC FIDFRs IREFLT»R END) 


INTEGRATES FROM REND (RAY BOTTOM) TO R3 
USING CENTER POIAT INTEGRATION FORMULA 


REAL*8 R(L$)yV(L$) 
REAL¥8 RADDEGsREsR1,R2-R3 sVREDe TINGS AIL sAI2Z,AINC DAI, ATs DELTA, TIME 
REAL*8 DELTAC,DELTAM,DELTAP)EMy EPs SLOPE, PyRAD yDELTAKyDELTAD,TRED 
REAL SA V1 }PEZ+TESTs DR, RADC,RADT, VELCs VELT»ROVsA,THOLD, VEL, REND 
EURMAUTN TD Ah od Kat XP IX" Ze 5 (Xe 818) ) 
FORMAT(8(1X»D14.8)) 
Ic =1 
1H=0 
OCALL U VEL(R VoL Sy 
LRAD»sVEL,IC,IDIM) 
LE (PL EQ. 0.000 )P=(RAD/VELI*DSIN(AT) 
PE2=Px*xP 
IF(R32.26Q.0.0N0IR3Z=R(1) 
i= 
FESTHANDELCATVITEST) 
DR=1L/DFLOATC(IOR) 
IFC IOUT-EQe1L)WRITE(6,100) 


MODIFY FIRST STEP IN INTEGRATICN TO THIN HOMOGENEOUS SHELL FORMULA 
USE THIS FORMULA ONLY IF IREFLT =0 
OTHERWISE RETAIN PRESENT FORM 


CONT INUE 
RADC=RADFDR/2~.0D0 
RADT=RAD—DR¥TEST 
OCALL U VEL(R Vol» 
IRADC,VELC,IC,IDI™) 
OCALL U VEL (ReVel $e 
LRADT,sVELT,»IT,»IDIM) 

Pete ADU/VEl Is ore Io) 10.1 
TPIRADT/ VELG.GT.PIGO TO 2 
[PUR AD+OR. GE~RSIGO 10.3 
ROV=RADC/VELC 
A=RNV*RDV—-PE2 
A=OR/DSQRT(RDV*ROV=PE2) 
DELTA=DELT A+P ¥A/RACLC 
TIME=TIME*ROV*A/VELC 
RAD=RAD+OR 

fet LOUIS EGS LV GALL QUT (DELTA, TIMEsR Ex TINCs THOLDyRADs 1010) 
GO TO 4 

CONTINUE 

CONTINUE 

[DR=IDR/2 
DR=1/DFLOATCICR) 

IT=1 

het) DRe GI won. 2G0. 10 4 
DR=1L/DFLOCAT( IORI) 
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TEST=RAD/OR 

GO FQ) % 

CONTINUE 

NR=P 3-RAN 

RADC=RAD+DR /2 .000 
OGAELCUTVER(IR Vial S+ 
LRANC»VELC»,IC,IDI™) 
RDV=RADC/VELC 

IF (RDV¥ROV=-PE2.LE.0.000)G0 TO 5 
A=DR/DSORT( ROV*RDV—-PE2 ) 
DELTA=DELTA+P *A/RATCC 
TIME=TIME+RDV*A/VELC 

RAD=RAD+DR 

CON DINGE 

PREPOUTZEOSECALL CUT (DELTAyTIME,REyTINC,s THCLDyRADy1Dy1) 
RETURN 
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SUBROUTINE DOWN TCR Vel $y 


LATOR YE,» RoyRay L0R1,1OR2)17EST,»IDIM, IOUT,» DELTA, TIMEyREyPsRAD,THOLDS 
2TINC -IOR,ITREFLT,REND) 


INTEGRATES FROM R1 TO REND (RAY BOTTOM) 
USING THIN HOMOGENEOUS SHELL FORMULA 


REAL*8 K(L$),V(L$) 

REAL*& V1,PE2,TEST,DRyRADC,RADT, VELCye VELT sRDV ey Ay THOLDy VEL 

REAL *8 DELTAC,DELTAM, DELTAP,EMyEPy SLOPEyP yPAD gDELT AkKy DEL TAD, TRED 
REAL*8 RADDEG,REaR1LeR2yaR3eVREDe TINCs ATL oAI2Z,AINC DAI yA »DELTAy TIME 
REAL*8 RHOLD»PR END,yDRI2,RADMORy RADPOR: PVC yPVCE 2 

hCG fa) 

LD = 6 

IDR2D02=1DR2/2 

DELTA = O.Of0 

TIME = 0.0D0 

be(Riwele. 06.000) RL = RC) 

RHOLD =R1 

CALL D VEL (RoVel$y,R1leVLi,ICeIDIM) 

Pei Poe EO. 0-000) P = RIZDSINVAII/ VI 

Rie aaa esc 

CALL BOTTOM (ReVyICIM,P,AT,REND,IREFLT) 

TES Pos DELOATAITTEST) 

IF(IREFLT .EQ. 1) TEST = C.OD0 


RAD = R1 
IDR = LOR) 
CONTINUE 


ORes= 1.0007 DE LCAT CI UR) 
NRD2 = DR/2.0D0 


I = 0 
CONTINUE 
Posse) 


RADG =" RAD =DRD2 

RADMDR =RHOLD — FLCAT(I)*DR 

IF(RADMDR -—TEST#OR .LE~. REND) GO TO % 
CALL D VEL (RyVylSeRADC,VELCy IC, IDIM) 
PVC = P * VELC 

PYCE? PVC * PVC 

DELTA DELTA+DARCOS (PVC/RAD)—DARCOS (PVC/RADMDR } 
TIME = TIME+( DSCRT (RADERAD-PVCE2) -DSGRT( RADMDR¥RADMOR=PVCE2))/VELC 
RAD = RADMDR 

GQ TO 1 

CONTINUE 

Te REPEL 2bOss 1) GOTO 2 

RHOLD = RAD 

It = 0 

IDR = IDR*¥2 

TeCLUR wb.) 1DR2D2)TES1=0-000 

Hotes «leew ADRe GO 1002 

CONTINUE 

DR =(RAD -RENDI/2.000 

DRO2=OR /2-0D0 

Shoe RACY = ORD? 
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CALL Dy VER (RyVelS,RADC »VELC.1IC,IDIM) 
PVC, = P¥ VELC 

PVGEZSPVN Ce ee PVC 
DELTA=DELTA#DARCOS(PVC/RAD) 
TIME=TIME*(DSQRT(RAD*RAD=PVCE2))/VELC 
LEER EC Ie bOs0dG0. TG. 
DELTA=DELTA-DARCOS (PVC/REND) 
TIME=TIME-DSORT (REND#REND=PVCE2)/VELC 
CONTINUE 

RAD = REND 

RETURN 

END 


» so 
fom 


_ 
rye 


, 


(Ay A) KAS 


Dee tee 


SUBROUTINE UP TCR V L$, 


PAT oR R2R2s TUR] LORZeITEST pLDI My LOUT, DEL TAyTIMNESREs? sRADs THOME 
ATINC »TOR,ITREFLT,REND) 


INTEGRATES FROM REND (RAY BOTTOM) TO R3 
USING THIN HOMOGENEOUS SHELL FORMULA 


REAL*8 R(LS),VIL$) 

REAL *8 RADDEG  »REgR1yR2gR3IeVREDeTINCs ATL sAL2Z,AINC sDATI SAT SOEL Tey LMS 
REAL *8 DEL TACs DELTAMy DELTAP, EMy EP y SLOPE oP RAD DELT AKy DELTADy TRED 
REAL*8 V1» PE29TESTyORy RADCyRADT, VELC» VELT sROVs Ay THOLD op VEL 

REAL *8 RHOLD,RENDs ORD2,RADMDR, RADPOR y PVC yPVCE2 

TEST = DFEOATCITEST) 


Loe= =) 
TOes &0 
ILDR=1DR2 


DR=1.-0D0/D0FLOAT(IOR) 

DRN2=0R/2.0D0 

RANC=RAD+DRD2 

CALLOUS VEL (R»V,L$sRADCeVELC, IC, IDIM) 
PeUPVG/RAD- LI. 120D0IGG TOe5 
RADPOR=RAD+4DR 

PVC=P*VELC 

PVCE2=PVC*¥PVC 

DELTA=DELTA+DARCOS (PVC/RADPDR ) 
TIME=TIME+(DSQRT(RADPDR*RADPDR-PVCE2) I/VELC 
RAD=RADPDR 

CONTINUE 

CALL °U VELIR GV el SsRADe VELv IC, 101M) 
Tele etes Os 0U0IGD 10 4 
P=RAD*DSIN(ALT)/VEL 

TREFLT=2 

CONT INUE 

TEC TREFLT © £Q-0)P=RAD/VEL 

LF(RAD/VEL eLT«P)PH=RAD/VEL 

PEZ2 = P*P 

(eset rea Oe ODO) Ra: = R(1) 

ie =O) 

CONTINUE 

RHOLD = RAD 

Nk = 1.O0DO/DFLOAT(IDR) 

HRo2 = DR/2s0D0 

CONTINUE 

RANC = RAD + DRD2 

Wa leat 

RADPDR =-RHOLD + FLOAT(I)*DR 

Te URADPOR «GE. R3) CONTE) 

uA eeaw VEL (Ry Veh $sRADCeVELCeICyIDIM) 
BYG. = Pe VELC 

PVCE2 = PVC * PVC 

DELTA =DELTA+DARCOS ( PVC/RADPOR)-DARCOS (PVC/RAD) 
DIMES TIMES {DS ORT ( RADPOR#RADPDR-PVCE2)-OSORT (RAD#RAD-PVCEZ} YN ELE 
PAD = RADPCR 

IF (RAC -TEST*OR GT. REND) GO TO z 
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PD Ress. LORE 


DR =e bs0D07 DF LGAT CIDR) 
Lest a=tRADLOR 

GaetO.3 

CONTINUE 

IOR = IORI] 

ORA= eh 49+) RAD 

DROZe= DR/2Z2.0D0 

RADC =(R34#RAD)/2.000 
CALL U VEL (ReVel$sRADC,VELC,IC,IOIM) 
PVG s="P* ©VELC 

PVGE2 PVC¥PVC 


DELTA = DELTAFt(DARCOS (PVC/R3)-—DARTOS (PVC/RAD) ) 
TIME = TIME#(DSQRT(R3¥*R3-PVCE2) -DSORT(RAD*¥RAD-PVCE2)) 


RAD = R3 
IDR = IDR 
RETURN 


END 
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OSUBPOUTINE UP I(ReVel$s 
LAL sRLeR2eR3e1TOR1L, IOR2 +1 TEST IDIMy IOUT ,DELTAsTIMEsREyP yRADs THOLD, 
ZV ING, LOR eI REPLT » REND 


INTEGRATES FROM REND (RAY BOTTOM) TO R3 

exe ASSUMING 83 = RI 

#KKASSUMING SAME VELCCITY STRUCTURF 4S USED IN DOWN I 

KR RENOT ERR RRR RE THIS ROUTINE DOFS NOT CALCULATE RAY PATH #*# 


REAL*8 R(L$),V(IL$) 
REAL*8 RADDEGsREyR1eR2yRBeVREDeTINCeAT1LyAI2,AINC DALAT DELTA, TIME 
REAL¥8 DELTACsDELTAM, DELTAP,EMyEPy SLCPE®P yRAD yDELT AKy DELTADy TRED 
REAL®B V1,PE2,TEST»DRyRACC »RACTyVELCy VELT sROVe Ay THOLDy VEL» REND 
100 FORMAT(1LXs507Xo*X* oe TX tZ 5 TXe*T8)) 
101 FORMAT(8(1X%,D14.8)) 
DELT A=DELTAFDELTA 
TIME=TIME+TIME 
RAD=R 1 
IDR=IDRI 
RETURN 
END 
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OSUBROUTINE D VEL(R yVyLS$y 
LRAD,VEL,IC ,IDIM) 


Ch Ch GY Oe) 


nN 


CALCULATES VELOCITY AT A GIVEN RADIUS FOR DECREASING R 
USING LINEAR INTERPOLATION FORMULA BETWEEN V( I), VOTt1) 
WHERE RI(I1).GT.R.GE-R(IT+1) 


REAL *8 R(L$),VIL$) 

REALS RADDEG RE aR] R2gR3eVREDeTINCsAL1eAI2,AINC yDAIy AI» DELTA, TIME 
Remus DEL TACs DELTAM,DELTAP, EM, EP, SLOPE, P yRADsDELT AKy DELTAD, TRED 
REAL*R Vl_PE2,TESTsDRyRADCyRALDTs VELCy VELT sROVe Av THOLD VEL 
IF(RAD.GE.R(1))9GO0 TC l 

ISTART=IC 

PDIML=TDIM=—1 

Pa TSISTARTs 1D IMI 

Ic=1 

LPS E+b 

PEXRAGVGE-R(IP1))GO TO 3 

CONTINUE 

CONTINUE 

VEL=CEV(TPLINV(T)) #RADAV( TD) *RCTPLIAVCIPLIFRITD7Z(RCIPLI~-R(10) 
RETURN 

CONTINUE 

VEL=V(01) 

Ic=1 

RETURN 

END 
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((1)a- (491 781) (1) 864 14) v= 81088 CVO aL EOI = | 


OSUPIRM UN EVES Ra gules 
TREAD) ee Viele eiCesal Del) 


GALOCULAT ES eVRLOGITY AT A, GIVEN RADIUS: FOR INCREASING BR 
USING LINEAR INTERPOLATION FORMULA RETWEEN VC), VOI+1) 
WHERE PUM) «Gti. GESRUI41) 


REAL*8 R(L$),V(L$) 

REAL*B RADDEG RE RL oK2yR3gVREDe TINC ALL xAI2eAINC DAI y AL yDEL TAs TIME 
REAL*B DELTAC, DELTAMy DELTAP,EMyEPy SLOPE®P RAD DELTAKy DEL TAD, TRED 
REAL*®8 V1,PEZ2/TESTsDR »RADC sRADT,VELC yVELT sROVeAgTHOLDy VEL 
TSTART=IC 

OO 4) T=ISTART,IDIM 

JiR Mae & 

[c=] 

PEW Ape eR CI Ge 10" 2 

CONTINUE 

CONTINUE 

TRGRAD @67.R(1))G0O 10 3 

T=J 

IPlL=J+1 

VEL= (UVC IPLINKV(T)) FRADAV (CL) #ROTPLIAVOTPLI*RITIIZORCIPLI-RULT)) 
RETURN 

CONTINUE 

VEL=V(1) 

RETURN 

END 
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100 


OSUBROUTING SET(RDEL sVDFL » JS eT Se RMAN,VMAN, KS yRo Vol Sy 


IMODELS 101) 
SETS ReV TO RDEL»VDEL sRMAN,VMAN ASSUMING NUMBER *®* MODEL 


REAL*8 RILS) VILE) 

RFAL*4 RDEL(J $5 1$),VOEL( J$,1 6) +.RMAN( KS) » VMAN( KS) 
RFAL*®R RADDEGyRE Rl yR2eR32eVREDATINC ALL eAILZeAINC pDAIL, AT, DELTA, TIME 
REALX® DELTAC sDELTAM,DELTAP, EMEP, SLOPEsPyRADy» DELT AK, DELTADyTRED 
REAL¥8 V1l_PE2eTEST »DRyRADC,RADT, VELCs VELT sROVeAyTHOLDS VEL 

Dork ® 1241 51S 

IF(RDEL(MODEL,1).LE.0.0)G0O TO 2 

[l=] 

R(LJ])=DFLOATCIFIX((RDEL(MCNEL,I11)*1.0E2)+0.5) )/1.0D2 

VI ITT) =DELOATC(IFIXC(VDELC MODEL, II) ¥*1.064)40.5))/1.0D4% 

CONTINUE 

CONTINUE 

Pl=1 T+! 

L$M1l=L$-1 

90 3 [T=I1,L$M1 

Aft aaa ee ol 

IF(RMAN(J).LE-0.0)G0 TC 4 

TOIM= 

R(1)=CFLOATCIFIX((RMAN(J)*1.20E2)+0.5)9/1.0D2 

V( IT) =DFELOATCIFIXCCVMAN(CS) ¥12-0E4)40.5))/1.004 

CONT INUE 

CONTINUE 

IDIM=IDIM+1 

R(1IDIM)=0.0D0 

VCIDIM)J=VCIOIM-1) 

FORMAT(iX,l0F13.5) 

RETURN 

END 
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CSUBROUTINE MOD(DELMOD, J$, 
LDELTAK MODEL) 


CALCULATES MODEL NUMBER (MODEL) FOR A GIVEN RANGE FROM SOURCE 


REAL*8 RADDEG,REsR1eR2eR3eVRECeTINC AIL yAI2,AINC yDATy AI, DEL TAs TIME 
REAL*®8 DEL TAC, DELTAM, DELT AP, EMyEPy SLOPE PeR ADs DELTAKs DEL TAD, TPED 
REAL*®B V1l,PFE2,TFESTyDRyRADCyRADT,VELCs VELT sROVe Ay THOLDYVEL 

RFAL*4 DELMOD(J$) 

DOL “Tt=1s5J3% 

MONEL=I 

FRUDEU TAK SE SDE MOD CI RETURN 

CONT TNUE 

RETURN 

END 
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SUBROUTINE OUT(DELTA,TIME,RE,TINCy THCLD,RAD,ID, I DUMP) 


DOU Toei Ay? PATH (CC-DORDINATES? X42 


MODY FOR CCALCDMP, PUGT! DUTPUT 


REAL ES 06015952 (5. T1590 

REAL*®8 RADDEG»sREsR1L eR2eR3eVREDeTINCAITLyAI2Z2,AINC »DAI,AI,DELTAg TIME 
REAL*8&8 DELTAC,DELTAM,DELTAP,EM,yEP, SLOPE, PyRADyDELT AKy DELTAD, TRED 
REAL *8 VL,PEZ,TESTsDRyRADC,RADT,VELC,VELT »RDV,A,THOLD VEL 
PURMATAAK. OU ZF belsFEs2)) 

LEGLDUMP weet JOU) FO 1 

IFCTIME=THOLD.LT-TINC)RETURN 

ID=iD+1 

X( ID =RAD* OS INI DELTA) 

Z(IDJ=RE-RADXDCOS( DELTA) 

TCIO) =TIME 

THO D=THOL DAT INC 

Pit Desi Day RETURN 

WRITE( 6,100) 00X0K) oe Z0(K),T(K)) sK=1,10) 

I N=0 

RETURN 

CONTINUE 

1D=1D+1 

X41 0) =RAD*¥DSINCODELTA) 

Z{ITD)Y=RE-RAD*XDCOS(DELTA) 

TED STIME 

WRITE( 6,100) 0(X(K)4Z¢(K),T(K) 9, K=1,ID) 

Tn=0 

RETURN 

END 
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Riese 
SUBROUTINE DATA(RDELyVDEL yDEL MOD, J $y 1 $¢RMANy VMANdKSoRo Vol $) 
RPEADSA VEUOCUTY STRUCTURE OATA CARDS SET It 


REAL¥*8 R(L$),VIL$) 

REALKG RDEL(J$,1$)¢VOEL(S $51 $) -DELMOD( I$) 9 RMAN(K$) 9 VMAN( KS) 

REAL¥BR RADDEG, RE gR1yR2yR3eVREDeTINCs ALL eAT2¢AINC DAI, AT DELTA, TIME 
REAL#R DELTAC,DELTAMyDELTAP,EMyEPy SLOPE yP RAD» DELT AKy DELTADy TRED 
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FORMAT(8F10.0) 
FORMAT (*0%,5(6Xe*RADIUS! e4Xe"VELNCITY')) 
FORMAT(1X,1L0F12.5) 
SNRMAT (*s MODEL'S 14," USED FOR DELTA LESS THAN ,FLOeLy *KM® ) 
EFORMAT(® MODEL® ,14,° USED FOR DELTA LESS THAN®,F14-8e°KM® ) 
FORMATCYILMODEL',15) 
WRITE(6,100) 
Oni id=1,50 


DATA CARDS SET 1-1 TERMINATETED WITH 9 IN COLUMN 80 
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CONTINUE 
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DAT e. CAR DS* SET® 242117 CARD)( £10.20 DELMOD(J) J=l INITIALLY 
AZIMUTH RANGE (KM) OF MODEL J LIES BETWEEN DELMOD(J-1),DELMOD(J) 


REAN(5,102)DELMOD(J) 
WRITE (6,105) J, DELMOD( J) 
DN 32 J=1,20%4 

N=1+3 


DATARICARDOS SET 1.3 8F10.0 ROEL,VDEL OF MODEL J 

TERMINATED BY READING 5 CARDS OR LAST RDEL OF CARD LE.O.-0 

TE RMAN (LAST OF CARD) TS CT 0.0 GO GN TO DATA CARDS “SET 16% 
TF RDEL CLAST OF CARD?) I5 0.0,READ NEXT DELMOM AND RMAN,VMAN 
BY READING NEXT SET OF DATA CARDS SET lels 1-2, AND 1.3 

THIS CAN BE DONF FOR 50 MODELS AT MOST 


READ (5,102) ((RDEL(JyK) pVDEL (J 9K)) 9K=T9N) 
IF (RDEL( JN) -LE20-.0)69 TO 4 

CONTINUE 

CONTINUE 

WRITE(69103) 

WRITE (65104) ((RDEL (JK) pVDEL( Jo K) Do K=19N) 
NMONDEL=J 

TE (RDEL(J,N)-LT-0-0160 TO 5 
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1 CONT ENUE 
ae CON TINGE 
DELMOD(NMODEL)=1.9E70 
WRITE (69106) NMQDEL, DELMOD(NMODEL) 


DATA CARDS SET 1-4 TERMINATED WITH 9 IN COLUMN 80 
VETIPE “OR MANTLE VEL CCl TY "STRUCTURE 


6 CONTINUE 
READL aS pLOvRrS Tee 
WRITE(6,101) 
IF(ISTOP.NE~9)GO TO 6 
NO 7 (7 98094 
N=[+3 


DATA CARDS SET 1.5 8F10.0 RMAN,VMAN 
TERMINATED BY READING 10 CARDS OR LAST RMAN OF A 


REAN (5,102) (( RMAN(K) » VMAN(K) ) 9 K=I9N) 
IF (RMAN(N).LE.0.90)G0 TO 8 
CONT INUE 
8 CONTINUE 
N=N=—-4 
N09 1[=1,54 
IN=I4N 
TF CRMAN(1I).LE.-0.0)G0 TO 10 
9 CONTINUE 
10 CONTINUE 
RMAN{IN)=0.0 
LE (VMAN( IN) «LE*000) VMAN(INJ=VMANCIN~1) 
WRITE(6,103) 
WRITE (69104) ((RMAN(K},VMAN(K) ) 9K=1lyIN) 
NO Li P=l,NMONEL 
OG ALE SETORCDEL yVOEL ydSe IT $e RMANy VMANo KS eRe Vol So 
LS Ce 8 
WRITE(6,107)1 
WRITE(6,103) 
WRITE (6,104) (CR(K) 2 VIK)) pK=l, TDIM) 
PL CONTINUE 
RETURN 
END 
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SUBROUTINE BOTTOM(R»V,IDIMyPyATsREND, IREFLT ) 


CALCULATES THE MINIMUM RADIUS (REND) GF A GIVEN RAY 


PREP = tOletNUGASE (OF REPRACTED RAY 
Pert sb tN GASE OF REELECTED RAY 


REAL*8 R(IDIM), VCICIM) 
REAL *8 PyAT,RENDsA185C 
[Pee Sy 293 

CONTINUE 

Ca" (ECL VFDSINTALINZV C1) 
CONTINUE 


TREELA = 0 
LODAMI = PDIM = 1 
BOT =" 14 1D 1M) 


Le aia 

feUe(TPl) ces. RU1)) GO VO 1 
B Ste wy BCL) 

A (VOIP) = VCI))/8 

B (Otte * VIII — RA LeV (TP 1078 
REND =(8*P)1/(1.0D0-(A*P)) 

[fe tREND «Ges RCLPL)) GO. TO: 5 
CONTINUE 

CONTINUE 

Ve CUREND! sll Rela) (GU 107.6 
IREFLT = 1 

REND = R(T) 

CONTINUE 

RETURN 

END 
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